# example.R
#
# example of how to work with the model and results from a wrapper language
# outside of GAMS. this example:
#    1. runs the baseline model
#    2. write a gams policy input file for the model that reduces the
#       productivity of capital in the electricity sector by 10%
#    3. reruns the models with the policy scenario
#    4. presents the change in key variables
#
# this example is intended to be run from the top level directory of the model
#

# load in the R utilities for sage
source(file.path("utilities","R_utilities.R"))



# ------------------------------------------------------------------------------
# run the model
# ------------------------------------------------------------------------------

# use the default aggregation benchmark file
benchmark_file = "default_aggregation.gdx"

# run the baseline
run.model(output_file="baseline.csv",benchmark_file=benchmark_file)

# create a gams input file containing the policy shock
policy_file = "example_shock.gms"
fc = file(file.path("examples",policy_file),open="wt")
writeLines("* new capital productivity shock in electricity sector",con=fc)
writeLines('prod_ind(t,r,"k","ele","new") = 1.1*prod_ind(t,r,"k","ele","new");',
           con=fc)
close(fc)

# run the policy scenario
run.model(gdx_baseline_file="baseline.gdx",output_file="policy.csv",
          policy_file=paste0("examples\\\\",policy_file),
          benchmark_file=benchmark_file)



# ------------------------------------------------------------------------------
# display the results
# ------------------------------------------------------------------------------

# List of parameters to output results for
vars = c("gdp","cd","ld","leisure")

# load the results
bau = read.csv(file.path("output","baseline.csv"),stringsAsFactors=FALSE)
policy = read.csv(file.path("output","policy.csv"),stringsAsFactors=FALSE)

# years in the analysis
years = unique(bau$time)
years = years[!is.na(years)]

# Run the comparison
for (v in vars) {
  for (t in years)
    compare(v,t)
}

# output equivalent variation
ev = sum(policy$value[policy$parameter=="ev"])
printf('EV........................: %7.3f',ev)
